// 
//  WgsSimulator/WgsSimulator.hh ---
// 
//  $Id: WgsSimulator.hh,v 1.6 2000/05/01 14:59:35 harley Exp $
// 
//
//  Author:  James H Gorrell <harley@bcm.tmc.edu>
//  Keywords: 
//
//  Commentary:
//  * James wrote the initital version.
//  * An Object to generate simulated shotgun reads.
//


// So we don't have to fuss with whether the file is included multiple 
// times.   I put this in all of my include files.   
#ifndef _WGSSIMULATOR
#define _WGSSIUMLATOR


#include "LocalSTL.hh"
#include "Sequence.hh"
#include "FASTAFile.hh"
#include "RandomNumberGenerator.hh"
#include "MiscUtils.hh"
#include "Sequence3Bit.hh"

class TraceQualityTemplate{
public:
  TraceQualityTemplate(){
	 mTqFileIdx = 0;
	 mTqIdx = 0;
	 mMinP20 = 20;
  }

  int mTqFileIdx;
  int mTqIdx;
	int mMinP20;

  vector<string> mTqFileList;
  SequenceQualityVector mTqVec;

  
  void GetNextTraceQuality(SequenceQuality &Tq);

	void GenerateTraceQuality(float errRate,SequenceQuality &Tq,int rsize);

  void ReadNames(string &FileName){
	 ReadFOF(FileName,mTqFileList);
  }
  
};


//---------------------------------------------------
// * WgsSimulator
//---------------------------------------------------
// 
class WgsSimulator{
public:

  WgsSimulator(){
	 mReadNum = 0;
	 mReadCounter=0;
	 mSeqIdx = 0;
	 mTotalSeqSize = 0;
	 mReadNameCount = 0;
	 mMinP20 = 20;
	 mErrRate = -1;
  }

	

  Sequence3BitVector mSourceSequences;

  RandomNumberGenerator mRNG;

  TraceQualityTemplate mTqTemplate;

  int mReadNameCount;

  int mSeqIdx;
  int mReadNum;
  int mReadCounter;
  long mStartTime;
  int mTotalSeqSize;

	int mMinP20;
	float mErrRate;
  
  
  void AddReads(string ClonePrefix,int NumClones, 
					 float MeanCloneSize, 
					 float StDevCloneSize,
					 bool bPaired,
					 ostream &fout,ostream &qout);
  
  void GenerateFwdReadWithErrors(int CloneStart,
											SequenceQuality &Tq,
											Sequence &Read,
											SequenceQuality &ReadQuality);

  void GenerateRevReadWithErrors(int CloneEnd,
											SequenceQuality &Tq,
											Sequence &Read,
											SequenceQuality &ReadQuality);

  char WgsSimulator::RandomBase();


  void SetReadName(Sequence &Read,string &prefix,string &suffix,char chem,int count);
  void SetReadName(SequenceQuality &ReadQuality,string &prefix,string &suffix,char chem,int count);


  /*
  // The read 'palette'...
  // ...Svec inserts
  int mCxtInsertLength;
  int mCxtInsertSdev;
  // ...Read params
  int mCxtReadAvgLength;
  int mCxtReadSdev;
  int mCxtReadMaxLength;

  // names
  string mReadNamePrefix;       // WGSAA
  int    mReadNameAttempt;      // 1
  int    mReadNameCnt;          // How many reads so far
  string mReadChem;             // dye, term

  float mCxtPartChimeric;
  float mCxtPartU;
  float mCxtPartURError;

  // Sequence 
  SequenceVector        mReads;
  SequenceQualityVector mReadsQual;

  // Stats
  int mTotalPhred20;
  int mCountReadUR;
  int mCountReadU;
  int mCountReadChimeric;
  int mCountReadURError;

  vector<float> mErrorRate;

  long mStartTime;

  // Statistics
  int mInsCount[3000];
  int mDelCount[3000];
  int mErrCount[3000];
  int mBaseCount[3000];
  
  // Options
  string   mPhrapCommand;
  int      mLongMean;
  float    mLongFraction;
  int      mReadStDev;
  float    mReverseFraction;

  // 
  RandomNumberGenerator mRNG;

  // tada!
  WgsSimulator() {
    // Default to what we generate
    mCxtInsertLength = 2000;
    mCxtInsertSdev   =  300;
    // ...Read params
    mCxtReadAvgLength   = 500;
    mCxtReadSdev        = 100;
    mCxtReadMaxLength   =   0; // Zero based

    // ...
    mReadNameCnt     = 1;
    mReadNamePrefix  = "WGSAA";
    mReadNameAttempt = 1;

    // perfect data
    mCxtPartChimeric=0;
    mCxtPartURError=0;
    // paired reads
    mCxtPartU=0;

    //
    mTotalPhred20   = 0;
    mCountReadUR    = 0;
    mCountReadU     = 0;
    mCountReadChimeric   = 0;
    mCountReadURError   = 0;

    // zero out stats 
    for(int i = 0;i<1500;i++){
      mInsCount[i]  = 0;
      mDelCount[i]  = 0;
      mErrCount[i]  = 0;
      mBaseCount[i] = 0;
    }
  };

  // *poof*
  ~WgsSimulator(){};

  // Methods

  // jhg
  int  SvecInsertSize();
  void GenInsertEnds(Sequence *, int &, int &);
  int  GenReadLength();
  void AddErrorsAndPushSeq(Sequence &);
  void GenaddReadUniversal(Sequence *);
  void GenaddReadUniversalReverse(Sequence *);
  void GenaddCntReadsUR(Sequence *, int, int);
  void GenaddCntReads(Sequence *, Sequence *, int, double);
  void printCxt();
  void printCounts();
  void GenaddReadChimeric(Sequence *, Sequence *);
  void SetReadName(string &,int,char);

  // james
  float InDelRate(int pos,int readLen);
  char RandomBaseError(char oldBase);
  char RandomBase();
  //int  PositionDependentPhredValue(int pos,int readLen);
  int  CrudePositionDependentPhredValue(int pos,int readLen);

  void AddSequencingErrors(Sequence &ErrorFreeSeq,
                           Sequence &SeqOut,
                           SequenceQuality &QualOut);
  void AddSequencingErrorsWithBounds(Sequence &, Sequence &, SequenceQuality &, float,float);

  void GetReadEndpoints(int ShortMean,long &start,long &end);
  void SimulateShotgunReads(float TargetCoverage,float TargetPhred20);
  void SimulateShotgunReads(int numReads,float TargetPhred20);
  void SimulateShotgunReads(float TargetCoverage,float TargetPhred20,
                            SequenceVector &Reads,
                            SequenceQualityVector &ReadsQual);

  void AddSimulatedShotgunReads(float TargetCoverage,float TargetPhred20,
                                SequenceVector &Reads,
                                SequenceQualityVector &ReadsQual);
  void SimulateShotgunReads(int numReads,float TargetPhred20,
                            SequenceVector &Reads,
                            SequenceQualityVector &ReadsQual);

  void SimulateShotgunReads2(float TargetCoverage,
                             float TargetPhred20,
                             SequenceVector &Reads,
                             SequenceQualityVector &ReadsQual,
                             vector<int> &StartLocs,
                             vector<int> &EndLocs);


  void SimulateErrorFreeShotgunReads(int numReads,int targetLength);
  void GenerateShotgunRead(int readNum,Sequence &Seq,int ShortMean);
  void GenerateShotgunRead2(int readNum,Sequence &Seq,int ShortMean,
                            long &ReadStart,long &ReadEnd);

  void RunSimulation();

  void SaveErrorStats(string FileName);

  void ReadPositionDependentErrorRates(string FileName);
  double PositionDependentErrorRate(int pos,int readLen);
 
  */

};

/*

//---------------------------------------------------
// * AsmSimStats
//---------------------------------------------------
// 
//
class JHGAsmSimStats{
public:
  
  JHGAsmSimStats() {
    mNumReads = 0;
    mReadSizeMean = 0;
    mPhred20Mean = 0;
    mPhred20Variance = 0;
    mContigs = 0;
    mContigSizeMax = 0;
    mActualCoverage = 0;
    mPhrapTime=0;
  };

  ~JHGAsmSimStats(){};

  // Simulation parameters
  vector<double> mErrRate;
  vector<double> mInsRate;
  vector<double> mDelRate;

  // Info about input
  int    mNumReads;

  int    mReadSizeMean;
  float  mPhred20Mean;
  float  mPhred20Variance;

  float  mActualCoverage;

  // Info about results. 
  int    mContigs;
  long   mContigSizeMax;

  int mPhrapTime;

  void PrintNumContigs(ostream &fout){
    fout<<mContigs;
  };

  void PrintContigSize(ostream &fout){
    fout<<mContigSizeMax;
  }
};


//------------------------------------------------------------
// * TableVec
//------------------------------------------------------------
// Row in an integer table
//
class JHGStatVec: public vector<JHGAsmSimStats>{
public:
  JHGStatVec(){};
  ~JHGStatVec(){};
};

//------------------------------------------------------------
// * StatTable
//------------------------------------------------------------
// A 2D integer table. (KJD: Of course, now I have the snazzy Table class
// which is more generic and feature rich).   
//
class JHGStatTable: public vector<JHGStatVec>{
public:
  JHGStatTable(){};
  ~JHGStatTable(){};
  
  void   Allocate(int rows,int cols);
  void   Clear();

  void   PrintNumContigs(int longMean,string FName);
  void   PrintContigSize(string FName);

  int    mRows;
  int    mCols;
};

*/ 

#endif
